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i Abstract 

The present work is concerned with the shape optimization of flapping wings in forward flight. The 
analysis is performed by combining a gradient-based optimizer with the unsteady vortex lattice 
method (UVLM). We describe the UVLM implementation and provide insights on how to select 
properly the mesh and time-step sizes to achieve invariant UVLM simulation results under further 
mesh refinement. Our objective is to identify a set of optimized shapes that maximize the propul- 
sive efficiency, defined as the ratio of the propulsive power over the aerodynamic power, under 

7-h lift, thrust, and area constraints. Several parameters affecting flight performance are investigated 

and their impact is described. These include the wing's aspect ratio, camber line, and curvature of 
the leading and trailing edges. This study provides guidance for shape design of engineered flying 

Z systems. 
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performance analysis, wing geometry description, camber effect, aspect ratio effect, leading and 
trailing edges curvature effects. 



1. Introduction 

The design and optimization of micro-air vehicles (MAVs) have been recently the topic of 
many research efforts |[I1^1H1SI|51IE1I71[S1I^1 HDl [IH H3 H31 US H5|. The MAV concept is 
quite new and the associated design space is large and not fully explored yet. These small flying 



Preprint submitted to Computers & Fluids 



November 13, 2012 



aeroelastic systems are expected to operate in urban environments and confined spaces (i.e., inside 
buildings, caves, tunnels) and endure a variety of missions such as inspection and observation of 
harsh environments inaccessible to other types of vehicles or where there is danger to human life. 
Consequently, their design must satisfy stringent performance requirements, such as high maneu- 
verability at low speeds, hovering capabilities, high lift to sustain flight, and structural strength 
to survive gust loads and possible impacts. These requirements can be achieved mainly through 
two propulsion mechanisms: rotating helicopter blades or flapping wings [IV]. Through observing 
and simulating insects and birds flights, it has been concluded that flapping wings offer greater 
efficiency, especially at small scales [fT8l [T9l l20ll . As a step toward designing efficient flapping- 
wing vehicles, several studies QEDIZ2I23I21|23I2S|23|23 explored and investigated the 
characteristics of these flying animals. 

Although it is well known that natural flyers exploit a variety of mechanisms and aerodynamic 
aspects to control and manoeuver their flights, experimental observations do not enable a complete 
understanding of the physical aspects and dynamics of flapping flight. As such, there is a need 
to model the unsteady aerodynamic aspects of flapping- wing vehicles. Computational modeling 
and simulation are necessary to evaluate the performance requirements associated with flapping 
flight and to identify the relative impact of different design parameters (e.g., flapping parameters, 
shape characteristics). Several computational modeling strategies that are based on variable fidelity 
physics have been reported in the literature IB 111 [2ll [30l SB E21 [33l [34l Ell . Willis et al. [EH devel- 
oped a multifidelity computational framework that involves a combination of potential flow models 
and a Navier-Stokes solver. They showed how the use of different levels of geometric and physical 
modeling fidelity can be exploited to ease the design process of flapping wing systems. Certainly, 
the higher-fidelity Navier-Stokes simulations incorporate a more complete physical model of the 
flapping flight problem, however, the extensive computational resources and time associated with 
the use of these tools limit the ability to perform optimization and sensitivity analyses in the early 
stages of MAV design. Thus, to alleviate this burden and enable rapid and reasonably accurate 
exploration of a large design space, it is fairly common to rely upon a moderate level of modeling 
fidelity to traverse the design space in an economical manner [@1[37||- As such, several research 
efforts have considered the use of the unsteady vortex lattice method (UVLM) for the design of 
avian-like flapping wings in forward flight [fTl[6ll7ll30ll38l[39ll40l. 

The present work is concerned with shape optimization of flapping wings in forward flight. 
Our approach to the problem combines a local gradient-based optimizer with UVLM. The opti- 
mizer is based on the globally convergent (to a stationary point, not necessarily a global solution) 
method of moving asymptotes (GCMMA) [|4Tll42l . It belongs to the class of sequential approxi- 
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mate optimization methods and employs conservative convex separable approximations for solving 
inequality constrained nonlinear programming problems. It searches by generating approximate 
subproblems at each iteration, in which both the objective and constraint functions are replaced by 
convex functions. The construction of these approximating functions is based mainly on gradient 
information at the current iteration point in the design space. 

The unsteady vortex lattice method computes the forces generated by pressure differences 
across the wing surface resulting from acceleration- and circulation-based phenomena. This ac- 
counts for unsteady effects such as added mass forces, the growth of bound circulation, and the 
wake. Since the formulation of UVLM requires that fluid leaves the wing smoothly at the trail- 
ing edge (through imposing the Kutta condition), it does not cover the cases of flow separation at 
the leading-edge and extreme situations where strong wing-wake interactions take place. In other 
words, UVLM applies only to ideal fluids, incompressible, inviscid, and irrotational flows where 
the separation lines are known a priori. Nevertheless, Persson et al. 113011 showed through a de- 
tailed comparison between UVLM and higher-fidelity computational fluid dynamics simulations 
for flapping flight that the UVLM schemes produce accurate results for attached flow cases and 
even remain trend-relevant in the presence of flow separation. As such, in [30J the authors rec- 
ommended the use of an aerodynamic model based on UVLM to perform preliminary standard 
design studies of flapping wing vehicles, specifically, for cruise configurations where there is no 
flow separation nor substantial wing-wake interactions that would degrade the performance of the 
vehicle. The reduced fidelity afforded by UVLM makes its associated simulation time for each 
flow configuration is on the order of minutes on a desktop. 

In this paper, we seek to identify a set of optimized shapes that maximize the propulsive ef- 
ficiency of a flapping wing under several constraints. We explicitly include lift, thrust, and area 
constraints. The wing geometry is described using B-spline parameterizations, which are the stan- 
dard technology for describing geometries in computer-aided design (CAD) Il43ll44ll45ll46l . This 
choice simplifies the design optimization process by enabling mesh generation directly from the 
CAD model. The associated basis functions can be used to smoothly discretize wing shapes with 
few degrees of freedom. The location of the control points constitutes the design variables. Results 
suggest that changing the shape yields significant improvement in the flapping wings performance. 
This study is our first step towards constructing a unified framework that will facilitate the design 
of engineered flying systems. 

The remainder of the paper is organized as follows: first, the aerodynamic model used to 
simulate flapping flights is detailed. The B-spline representation that defines the wing geometry 
is then introduced, followed by the formulation of the shape optimization problem. The optimal 
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shapes obtained for different flight configurations are reported and insights of the trailing and 
leading edges curvatures, camber, and aspect ratio effects on the flapping flight performance are 
presented. The work concludes with a discussion of the main features of the optimal shapes. 

2. Aerodynamic Modeling of Flapping Wings 

We use a three-dimensional version of the unsteady vortex lattice method (3D UVLM) to sim- 
ulate the aerodynamic response of flapping wings in forward flight. This aerodynamic tool is 
capable of simulating incompressible and inviscid flow past moving thin wings and capturing the 
unsteady effects of the wake, but not the viscous effects, flow separation at the leading-edge, and 
extreme situations with strong wing- wake interactions. Unlike standard computational fluid dy- 
namics schemes, this method requires meshing of the wing planform only and not of the whole 
flow domain. UVLM uses a global vorticity circulation balance, which is described on a vortex 
lattice. Given the circulation on a ring, UVLM relies on the Biot-Savart law to construct the ve- 
locity field at every point in 3D space. Using the wing surface discretization, UVLM imposes 
a collocated no-penetration condition to construct a potential flow about the mid- surface of the 
wing. The satisfaction of the no-penetration condition allows us to compute the circulation in each 
vortex loop, closing the system. By this sequence of approximations, UVLM reduces the demand 
for computational resources. These features make the method competitive for performing opti- 
mization studies that require many simulations by appropriately capturing flow dynamics at a low 
computational cost. 

Features of the UVLM solver used in this study include 

• The shape of the wing is generated using a B-spline representation. 

• The wing surface is discretized into a lattice of vortex rings. Each vortex ring consists of 
four short straight vortex segments, with a collocation point placed at its center. 

• For rigid wings, the position of the grid points is specified by applying a sequence of rotations 
(pitching and flapping). 

• The no-penetration condition is imposed at collocation points. The normal component of the 
velocity due to wing-wing interactions, wake-wing interactions, free-stream velocities, and 
wing rotations is assumed to vanish at each collocation point. Using the Biot-Savart law to 
compute velocities in terms vorticity circulations T yields a linear system of equations that 
can be expressed as: 

A wi ~ wi r wi = -A wa ~ wi r wa + v n (i) 
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where ^4 W1 ~ W1 and yt wa ~ wl are wing-wing and wake-wing influence matrices, respectively. 
The vector V n collects the normal component of the velocity at each collocation point due to 
the wing motion. The vectors T wi and T wa stand for the circulations of the vortex elements 
on the wing and wake, respectively 03 |6]]. The LAPACK library BTTl is used to solve the 
linear system given by Equation ([T]). 

• Vorticity is introduced to the wake by shedding vortex segments from the trailing edge. These 
vortices are moved with the fluid particle velocity and their individual circulation remains 
constant (i.e., T^ At = T^ a ). The wake elements have been truncated in the flowfield and 
load computation and only those which are located within 10 c are accounted for, where c 
is the chord length. This truncation was observed to speed up significantly the simulation 
while introducing negligible loss in the solution accuracy. 

• The pressure is evaluated at each collocation point based on the unsteady Bernoulli equation 
and then integrated over the wing surface to compute the aerodynamic forces and power. 

Further details of the derivation, implementation, and verification of this method and aerody- 
namic loads computation are provided in [Til[6l l30l[35ll48ll49l . 

3. Modeling of Wing Shape: B-spline Representation 

3.1. B-spline curves and surfaces 

B-splines are piecewise polynomial functions generalizing the Bernstein-Bezier polynomials. 
The B-spline basis functions of degree p, denoted by N i>p (^), associated to a non-decreasing set of 
coordinates, called the knot vector X — {£i, £ 2 , • • • , £n+ P +i}, are defined recursively as 



1, 

0, otherwise 



for i — 1 , . . . , n and p > 1. Knot multiplicities reduce the continuity of the basis at the location of 
the multiplicity. If a knot multiplicity of k is used, the continuity on the basis is C p ~ k at that knot. 
The basis becomes interpolatory at knots with multiplicity p whereas knot multiplicity of p + 1 
makes the basis discontinuous. 
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The B-spline curve of degree p with control points Pi, ... , P n is defined on the interval [a, b] = 
[£ p +i, Cn+i] as the linear combination of the control points and basis functions, that is, 

n 

i=i 

The piecewise linear interpolation of the control points is called the control polygon. A feature 
of B-spline curves is that the curve defined by the basis and control points lies inside the convex 
hull of the control polygon. This makes the control polygon useful for approximating the rough 
character of a curve. 

A B-spline surface is defined using tensor products of B-spline basis functions written in two 
parametric coordinates £ and 77. Let N i)P and Mj >q denote basis functions of degree p and q as- 
sociated to the knot vectors X = {£1, £ 2 , • • • , £n+ P +i} and y = {rj U r) 2 , . . . , r] m+q+1 }, and P ijt 
i = 1, . . . , n, j = 1, . . . , m is a net of control points in three-dimensional space. A B-spline 
surface is then defined as 

n m 
i=l j=l 

B-splines have long been used in the computer aided design community to model curves and 
surfaces. The reader is referred to l|43l 1441 |45ll46l for more details on B-splines. 

3.2. Wing shape parametrization 

In Figure [TJ we plot the wing geometry model based on a B-spline representation. The control 
points are shown as a network of connected spheres and the wing as a shaded surface. In B-splines, 
the surface does not interpolate the control points, rather the control polygon is an approximation 
to the surface which it discretizes. Perturbing the control points changes the shape of the wing. 
However, we deform the wing shape in such a way that it preserves some of its original features. 
In particular, we allow the camber edge to scale and translate in the y-direction only, as indicated 
by the blue spheres and arrows in the figure. The green spheres and arrows are similar except that 
we also allow translation in the ^-direction. The black sphere is held fixed in all cases. Finally, we 
proportionally redistribute the nodes to preserve to some extent the aspect ratio of each element of 
the mesh. We deform the wing surface in this way for three reasons. 

1. The camber edge is typically something known from experience and we do not want to 
depart completely from established designs. 
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2. Changing the camber significantly may lead to cases where, physically, flow separation oc- 
curs. This is a case we wish to avoid in practice as well as in our numerical method. 



The approach reduces the number of degrees of freedom such that the optimization is more 
computationally efficient. While Figure[T]uses 48 variables (16 control points in 3D space), 
we define our deformed wing using only 8 degrees of freedom (far end scales only, the 
middle rows scale and translate, and the near end scales, and translates in two directions). 



As an example, Figures |l(b) and |l(d) depict a possible deformed wing geometry based on our 
rules. We modify the wing shape in this way in all simulation results reported in this work unless 
stated otherwise. 




(a) Undeformed shape (b) Deformed shape 




(c) Undeformed shape - top view (d) Deformed shape - top view 



Figure 1: Wing geometry model based on B-spline representation. Spheres are used to represent the control points and 
arrows denote their perturbation directions. The straight line going over the spheres on the grid is the control polygon 
while the shaded surface represent the actual wing surface. 
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4. Shape Optimization 

4.1. Baseline wing shape 

In this section, we consider cambered rectangular wings with an aspect ratio of six (baseline 
case). The cambered wing has a NACA 83XX cross-sectional profile as studied previously by 
Stanford and Beran [1J and later by Ghommem et al. (6). The indication of the last two digits 
of the NACA profile, describing the wing thickness, is irrelevant for our study since the current 
implementation of UVLM assumes a thin wing. The symmetric flapping motion (about the wing 
root) is prescribed as: 

4>(t) = A^cos^u t), (2) 

where is the flapping angle and the flapping amplitude is set equal to 45°. Furthermore, the 
wing root is placed at a fixed angle of attack (pitch) of 5°. Different reduced frequencies are used 
to examine the effect of flow unsteadiness on the aerodynamic performance of flapping wings. The 
reduced frequency is defined as 

u c 

where tu is the flapping frequency, c is the chord length, and is the freestream velocity. The 
transient variations of the lift and thrust over one flapping cycle predicted by the current UVLM 
and those obtained by Stanford and Beran [[TJ| are shown in Figure |2j A good agreement can be 
clearly observed. For the sake of comparison, six elements are used along the chordwise direction 
and 10 are used along the semispan, and the flapping period is discretized into 40 time steps. In 
the next section, we show the need to consider smaller time steps and more refined aerodynamic 
meshes to obtain convergence. Our analysis provides guidance to the mesh size and time step size 
that need to be used to achieve convergent UVLM simulations under refinement. To our knowl- 
edge, the analysis presented in the next section is the first that achieves solution independence to 
discretization for UVLM. The rectangular flapping wing case, as described above, constitutes the 
baseline case of the optimization studies conducted next. 

4.2. Convergence analysis 

We first investigate the effect of mesh refinement on the lift coefficient while keeping the num- 
ber of time steps per flapping cycle N ts constant and equal to 40. Figure [3] shows that UVLM fails 
to converge as the mesh is refined. 

To fix the convergence issue, we perform several simulations while varying the number of time 
steps per flapping cycle N ts over a range of 40 to 160 and compare the transient variations of the 
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Stanford and Beran 201 
Current UVLM 




(a) Lift coefficient 



Stanford and Beran 2010 
Current UVLM 
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t/T 

(b) Thrust coefficient 



0.8 



Figure 2: Lift and thrust computed from UVLM for one flapping cycle: comparison with results obtained by Stanford 
and Beran [ 1 1. Six elements are used along the chordwise direction and ten are used along the spanwise direction (for 
half wing). 
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Figure 3: Transient variations of lift coefficient for different aerodynamic mesh sizes (fixed number of time steps per 
flapping cycle, N ts = 40). 

lift coefficient to select properly the time step size. The obtained results for a coarse mesh of six 
chordwise elements and 10 spanwise elements are shown in Figure |4j UVLM converges with a 
number of time steps per flapping cycle of 120. There is little change between the simulation using 
a number N ts of 120 and that using 160 time steps. 





A, N ts = 40 








V N ,s = 80 








\v / N =120 











0.2 0.4 0.6 0.8 1 

t/T 



Figure 4: Transient variations of lift coefficient for different time step sizes (fixed spatial resolution, 6x 10 mesh). 

As for the effect of the aerodynamic mesh on the prediction of the aerodynamic loads, the 
transient variations of the lift coefficient are evaluated for different chordwise and spanwise re- 
finements while keeping N ts constant and equal to 120. Results of the lift coefficient obtained for 
different numbers of chordwise and spanwise elements are shown in Figure [5] Clearly, refining the 
mesh along the chordwise direction affects substantially the prediction of the aerodynamic loads. 
We observe that the use of 24 chordwise panels yields convergence of the lift coefficient. The 
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0.4 0.6 

t/T 



0.8 



0.2 



0.4 0.6 

t/T 
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(a) Chordwise refinement 



(b) Spanwise refinement 



Figure 5: Transient variations of lift coefficient for different chordwise and spanwise refinements of the aerodynamic 
mesh. 

effect of the spanwise refinement on the convergence trend of the lift force is much less significant, 
as shown in Figure 5(b)| 

We propose to first determine the appropriate time step in order to observe convergence of 
UVLM simulations when refining the aerodynamic mesh. For the reminder of this work, the flap- 
ping cycle is discretized into 120 time steps and an aerodynamic mesh of 24 chordwise elements 
and 20 spanwise elements is used. We note that the use of 120 time steps per flapping cycle 
was checked for the aforementioned refined mesh size and was found appropriate to reach solution 
convergence; that is, as N ts is increased, the simulation results remain unaffected (not shown here). 

4.3. Problem formulation 

Wing shape selection presents a critical determinant of performance in flapping flight. In this 
study, we conduct series of shape optimization studies aiming at identifying a suitable set of wing 
planforms that enables efficient flights. The optimization problem is formulated as follows: 

Maximize r)(x) 
Subject to L*(x) > L* hh 

T*(x) > T* u , 

A(x) < A hh 



(3) 



max \ 9ij(x) 

y 

x e x, 



2 I — "cr? 



where r](x) is the propulsive efficiency, defined as the ratio of the propulsive power over the aero- 
dynamic power 03 HI, x is the vector of design parameters describing the perturbations introduced 
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to the locations of the control points, 



L* 



2L 



IT 



and P* 



IP 



are the normalized lift, thrust, and aerodynamic power, respectively, A is the wing's area, p is the 
fluid density, Uoo is the freestream velocity, the subscript 'bl' refers to the baseline case, and the 
overline denotes a time-averaged quantity over a flapping cycle. The wing geometry is constrained 
by restricting the feasible perturbations of the locations of the control points to lie in a closed and 
bounded set X. We chose to bound each component in x to vary in the interval [— c/2, c/2]. To 
illustrate the wing geometry parametrization based on B-splines, we plot in Figure [6] the aerody- 
namic mesh (6 x 10 lattice) and the location of the control points, represented by red circles, used to 
define the wing geometry. As pointed out in the previous section, the surface does not interpolate 
the control points and the aerodynamic mesh does not coincide with the control polygon. Here, 
the main objective is to maximize the cycle- averaged propulsive efficiency of the wing under lift 
and thrust constraints. The area of the wing is also restricted to be smaller or equal to a baseline 
value (i.e., the area of a cambered rectangular wing with an aspect ratio of six). In addition, we 
require that the angles 0y of each element i where j = {1, 2, 3, 4} that correspond to four edges 
of the element, in the aerodynamic mesh, see Figure [6| do not deviate too much from right angles. 
Large distortions of the grid could degrade the predictive capability of the aerodynamic model 
and also may give rise to undesirable unsteady flow effects such as flow separation and substantial 
wing-wake interactions. In our simulations, we specify the parameter 9 cr controlling the level of 
distortion as 9 cr = 15°. 

The optimization problem is solved with the globally convergent method of moving asymptotes 
(GCMMA) [|4Tll42l . The algorithm is supplied with numerically-computed gradients based on the 
first-order backward Euler scheme for both the objective function and the constraints. Several 
optimization runs are conducted where we vary the number of design variables and the polynomial 
degree of the basis functions employed in the B-spline representation. These simulations elucidate 
the effect of the design parameters on the numerical results. 

4.4. Results and discussion 
4.4.1. Single knot span case 

In the subsequent analysis, we use an aerodynamic mesh of 24 chordwise elements and 20 
spanwise elements and the flapping cycle is discretized into 120 time steps. We first consider the 
case where the B-spline representation is based on a single knot span [|43l|44l|45l|46]|. Thus, the 
number of control points increases with the polynomial degree. The optimal design from one run 
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Figure 6: Aerodynamic grid of the wing modeled using B-spline. Red circles represent the locations of the control 
points used to define the wing geometry. 

is used as the initial guess for the next run with a higher polynomial degree. The shapes obtained 
from the sweep of the optimization studies for different polynomial degrees are shown in Figure[7] 
The control points describing the optimal shapes are listed in Tables [l]-[4j Since the quartic and 
quintic responses were found similar, the geometric data is not reported for the latter. Optimal 
results show similar shape trends. For all cases, we observe that less area is distributed to the outer 
part of the wing. 

Table [5] provides a summary of the efficiency 77, lift L*, thrust T*, aerodynamic power P*, and 
area ratio A/Au obtained for the optimal configurations for different polynomial degrees. The 

Table 1 : Control points for the optimal shape: single knot span case and linear approximation is used in the B-splines 
representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


J 


X 


V 


z 


* 3 


X 


y 


z 








-0.015268 


0.000000 


0.000000 


1 


0.293877 


3.500000 


0.000000 





1 


0.153943 


0.000000 


0.089250 


1 1 


0.409020 


3.500000 


0.089250 





2 


0.323154 


0.000000 


0.126467 


1 2 


0.524162 


3.500000 


0.126467 





3 


0.492366 


0.000000 


0.034380 


1 3 


0.639304 


3.500000 


0.034380 





4 


0.661577 


0.000000 


0.091453 


1 4 


0.754447 


3.500000 


0.091453 





5 


0.830789 


0.000000 


0.031200 


1 5 


0.869589 


3.500000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


1 6 


0.984732 


3.500000 


0.000000 
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Table 2: Control points for the optimal shape: single knot span case and quadratic approximation is used in the 
B-splines representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


3 


X 


y 


z 


i 


3 


X 


y 


z 








-0.238849 


0.000000 


0.000000 


2 





0.740752 


3.500000 


0.000000 





1 


-0.032374 


0.000000 


0.089250 


2 


1 


0.834122 


3.500000 


0.089250 





2 


0.174101 


0.000000 


0.126467 


2 


2 


0.927493 


3.500000 


0.126467 





3 


0.380576 


0.000000 


0.034380 


2 


3 


1.020860 


3.500000 


0.034380 





4 


0.587050 


0.000000 


0.091453 


2 


4 


1.114230 


3.500000 


0.091453 





5 


0.793525 


0.000000 


0.031200 


2 


5 


1.207600 


3.500000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


2 


6 


1.300970 


3.500000 


0.000000 


1 





0.176031 


1.500000 


0.000000 












1 


1 


0.301366 


1.500000 


0.089250 












1 


2 


0.426702 


1.500000 


0.126467 












1 


3 


0.552037 


1.500000 


0.034380 












1 


4 


0.677372 


1.500000 


0.091453 












1 


5 


0.802708 


1.500000 


0.031200 












1 


6 


0.928043 


1.500000 


0.000000 













Table 3: Control points for the optimal shape: single knot span case and cubic approximation is used in the B-splines 
representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


3 


X 


y 


z 


i 


3 


X 


y 


z 








0.004282 


0.000000 


0.000000 


2 





0.112514 


2.000000 


0.000000 





1 


0.170235 


0.000000 


0.089250 


2 


1 


0.219737 


2.000000 


0.089250 





2 


0.336188 


0.000000 


0.126467 


2 


2 


0.326960 


2.000000 


0.126467 





3 


0.502141 


0.000000 


0.034380 


2 


3 


0.434183 


2.000000 


0.034380 





4 


0.668094 


0.000000 


0.091453 


2 


4 


0.541405 


2.000000 


0.091453 





5 


0.834047 


0.000000 


0.031200 


2 


5 


0.648628 


2.000000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


2 


6 


0.755851 


2.000000 


0.000000 


1 





-0.315068 


1.000000 


0.000000 


3 





0.491455 


3.500000 


0.000000 


1 


1 


-0.119715 


1.000000 


0.089250 


3 


1 


0.600877 


3.500000 


0.089250 


1 


2 


0.075637 


1.000000 


0.126467 


3 


2 


0.710300 


3.500000 


0.126467 


1 


3 


0.270989 


1.000000 


0.034380 


3 


3 


0.819722 


3.500000 


0.034380 


1 


4 


0.466341 


1.000000 


0.091453 


3 


4 


0.929145 


3.500000 


0.091453 


1 


5 


0.661693 


1.000000 


0.031200 


3 


5 


1.038570 


3.500000 


0.031200 


1 


6 


0.857045 


1.000000 


0.000000 


3 


6 


1.147990 


3.500000 


0.000000 
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Table 4: Control points for the optimal shape: single knot span case and quartic approximation is used in the B-splines 
representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


J 


X 


y 


z 


i 


3 


X 


y 


z 








0.089656 


0.000000 


0.000000 


3 





0.149637 


2.250000 


0.000000 





1 


0.241380 


0.000000 


0.089250 


3 


1 


0.249139 


2.250000 


0.089250 





2 


0.393104 


0.000000 


0.126467 


3 


2 


0.348641 


2.250000 


0.126467 





3 


0.544828 


0.000000 


0.034380 


3 


3 


0.448142 


2.250000 


0.034380 





4 


0.696552 


0.000000 


0.091453 


3 


4 


0.547644 


2.250000 


0.091453 





5 


0.848276 


0.000000 


0.031200 


3 


5 


0.647145 


2.250000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.746647 


2.250000 


0.000000 


1 





0.026628 


0.750000 


0.000000 


4 





0.444682 


3.500000 


0.000000 


1 


1 


0.165262 


0.750000 


0.089250 


4 


1 


0.555822 


3.500000 


0.089250 


1 


2 


0.303896 


0.750000 


0.126467 


4 


2 


0.666962 


3.500000 


0.126467 


1 


3 


0.442530 


0.750000 


0.034380 


4 


3 


0.778102 


3.500000 


0.034380 


1 


4 


0.581163 


0.750000 


0.091453 


4 


4 


0.889242 


3.500000 


0.091453 


1 


5 


0.719797 


0.750000 


0.031200 


4 


5 


1.000380 


3.500000 


0.031200 


1 


6 


0.858431 


0.750000 


0.000000 


4 


6 


1.111520 


3.500000 


0.000000 


2 





-0.292502 


1.500000 


0.000000 












2 


1 


-0.071887 


1.500000 


0.089250 












2 


2 


0.148727 


1.500000 


0.126467 












2 


3 


0.369341 


1.500000 


0.034380 












2 


4 


0.589956 


1.500000 


0.091453 












2 


5 


0.810570 


1.500000 


0.031200 












2 


6 


1.031180 


1.500000 


0.000000 
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Baseline Linear Quadratic Cubic Quartic Quantic 

Figure 7: Optimal wing shapes (single knot span case). Coordinates of the control points parameterizing the optimal 
shapes are reported in Tables [TJ|4] 

number of design variables is denoted by N D y. The optimal results show that changing the wing 
shape allows to achieve higher lift and propulsive efficiency but also more power is needed to 
be introduced to the flying system. As expected, rising the polynomial degree (i.e., increasing 
the number of degrees of freedom of the parametrization N^y) yields optimal shapes that enable 
flapping flights with higher efficiencies. Introducing swept and taper to flapping wings, which is 
obtained by using of linear polynomials in the B-spline representation, ameliorates significantly 
the flight performance. The use of higher polynomial degrees enables further improvements in 
the propulsive efficiency which stabilizes at fourth degree. Furthermore, we observe that the area 
constraint is active for all cases. Thus, there is a natural trade-off between the wing area and 
efficiency, but not so for thrust and efficiency. This indicates that, in the optimized configuration 
obtained at the low reduced frequency (k = 0.1), changing the wing shape enables to improve 
thrust generation without affecting the wing area. 

4.4.2. Fixed number of degrees of freedom case 

In this case, we keep the number of design variables (degrees of freedom) fixed and equal to 
iV D v = 12 and run the optimizer while varying the polynomial degree using in all instances the 
baseline configuration as the optimization starting point. The choice of N^y is made so that it coin- 
cides with the number of design variables required for the single knot span case when considering 
a fifth-degree polynomial. The corresponding optimal shapes are shown in Figure [8] The control 
points describing the optimal shapes are listed in Tables |6]-[8j Again, since the cubic and quartic 
responses were found similar, the geometric data is not reported for the latter. Figure [8] shows that 
all optimal shapes are qualitatively similar. In particular, the ones obtained using higher-degree 
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Table 5: Baseline vs. optimal results obtained for the single knot span case (k = 0.1). Lift, thrust, and area constraints 
are imposed. 



Wing shape 






V 


17 




P* A/A u 


Baseline shape 







0.219 


4.484 


0.185 


0.844 1 


Optimal 


Linear 


4 


0.426 


5.421 


0.4451 


1.0435 1 




Quadratic 


6 


0.433 


5.448 


0.415 


0.9566 1 




Cubic 


8 


0.439 


5.489 


0.443 


1.009 1 




Quartic 


10 


0.440 


5.470 


0.459 


1.040 1 




Quintic 


12 


0.440 


5.470 


0.459 


1.040 1 



polynomials look smoother. This smoothness in the shape affects the performance of the flapping 
wing. Although the single knot span case using fifth-degree polynomial and the fixed number of 
degrees of freedom cases have the same number of design variables, different shapes have been 
obtained. This difference can be associated to the choice of the starting shape when running the 
optimizer. Furthermore, we observe that the curvatures of the leading and trailing edges are both 
tilted up near the wing tip. Several studies have pointed out the importance of the curvature of the 
leading edge as being the location of the first contact of the fluid with the body and the camber 
line EES. 

Flight direction 

((T(( (( 
(( (( (( (( 

Baseline Linear Quadratic Cubic Quartic 

Figure 8: Optimal wing shapes (low reduced frequency, k = 0.1, and fixed number of degrees of freedom case, 
^Vdv = 12 )■ Coordinates of the control points parameterizing the optimal shapes are reported in Tables[6}j8] 

Table [^provides a summary of the efficiency rj, lift L*, thrust T*, aerodynamic power P*, and 
area ratio A/Au obtained for the optimal configurations for different polynomials. Clearly, linear 
polynomials do a good job in improving the performance of flapping wings. The optimal results 
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Table 6: Control points for the optimal shape: fixed number of degrees of freedom case and linear approximation is 
used in the B-splines representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


3 


X 


V 


z 


i 


J 


X 


y 


z 








0.019323 


0.000000 


0.000000 


3 





-0.160116 


1.800000 


0.000000 





1 


0.182769 


0.000000 


0.089250 


3 


1 


-0.009644 


1.800000 


0.089250 





2 


0.346216 


0.000000 


0.126467 


3 


2 


0.140827 


1.800000 


0.126467 





3 


0.509662 


0.000000 


0.034380 


3 


3 


0.291299 


1.800000 


0.034380 





4 


0.673108 


0.000000 


0.091453 


3 


4 


0.441771 


1.800000 


0.091453 





5 


0.836554 


0.000000 


0.031200 


3 


5 


0.592243 


1.800000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.742714 


1.800000 


0.000000 


1 





0.136362 


0.600000 


0.000000 


4 





-0.040048 


2.400000 


0.000000 


1 


1 


0.260379 


0.600000 


0.089250 


4 


1 


0.108170 


2.400000 


0.089250 


1 


2 


0.384396 


0.600000 


0.126467 


4 


2 


0.256389 


2.400000 


0.126467 


1 


3 


0.508413 


0.600000 


0.034380 


4 


3 


0.404607 


2.400000 


0.034380 


1 


4 


0.632430 


0.600000 


0.091453 


4 


4 


0.552825 


2.400000 


0.091453 


1 


5 


0.756447 


0.600000 


0.031200 


4 


5 


0.701043 


2.400000 


0.031200 


1 


6 


0.880463 


0.600000 


0.000000 


4 


6 


0.849262 


2.400000 


0.000000 


2 





-0.015657 


1.200000 


0.000000 


5 





0.244208 


3.500000 


0.000000 


2 


1 


0.112085 


1.200000 


0.089250 


5 


1 


0.389343 


3.500000 


0.089250 


2 


2 


0.239827 


1.200000 


0.126467 


5 


2 


0.534479 


3.500000 


0.126467 


2 


3 


0.367569 


1.200000 


0.034380 


5 


3 


0.679614 


3.500000 


0.034380 


2 


4 


0.495311 


1.200000 


0.091453 


5 


4 


0.824750 


3.500000 


0.091453 


2 


5 


0.623053 


1.200000 


0.031200 


5 


5 


0.969885 


3.500000 


0.031200 


2 


6 


0.750795 


1.200000 


0.000000 


5 


6 


1.115020 


3.500000 


0.000000 



18 



Table 7: Control points for the optimal shape: fixed number of degrees of freedom case and quadratic approximation 
is used in the B-splines representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


3 


X 


y 


z 


i 


J 


X 


y 


z 








0.060584 


0.000000 


0.000000 


3 





-0.223082 


1.875000 


0.000000 





1 


0.217153 


0.000000 


0.089250 


3 


1 


-0.071007 


1.875000 


0.089250 





2 


0.373723 


0.000000 


0.126467 


3 


2 


0.081068 


1.875000 


0.126467 





3 


0.530292 


0.000000 


0.034380 


3 


3 


0.233144 


1.875000 


0.034380 





4 


0.686861 


0.000000 


0.091453 


3 


4 


0.385219 


1.875000 


0.091453 





5 


0.843431 


0.000000 


0.031200 


3 


5 


0.537294 


1.875000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.689369 


1.875000 


0.000000 


1 





0.103366 


0.375000 


0.000000 


4 





-0.034651 


2.625000 


0.000000 


1 


1 


0.236516 


0.375000 


0.089250 


4 


1 


0.119043 


2.625000 


0.089250 


1 


2 


0.369666 


0.375000 


0.126467 


4 


2 


0.272736 


2.625000 


0.126467 


1 


3 


0.502816 


0.375000 


0.034380 


4 


3 


0.426430 


2.625000 


0.034380 


1 


4 


0.635966 


0.375000 


0.091453 


4 


4 


0.580123 


2.625000 


0.091453 


1 


5 


0.769116 


0.375000 


0.031200 


4 


5 


0.733817 


2.625000 


0.031200 


1 


6 


0.902266 


0.375000 


0.000000 


4 


6 


0.887511 


2.625000 


0.000000 


2 





0.031374 


1.125000 


0.000000 


5 





0.191797 


3.500000 


0.000000 


2 


1 


0.145334 


1.125000 


0.089250 


5 


1 


0.344134 


3.500000 


0.089250 


2 


2 


0.259295 


1.125000 


0.126467 


5 


2 


0.496471 


3.500000 


0.126467 


2 


3 


0.373256 


1.125000 


0.034380 


5 


3 


0.648809 


3.500000 


0.034380 


2 


4 


0.487216 


1.125000 


0.091453 


5 


4 


0.801146 


3.500000 


0.091453 


2 


5 


0.601177 


1.125000 


0.031200 


5 


5 


0.953483 


3.500000 


0.031200 


2 


6 


0.715137 


1.125000 


0.000000 


5 


6 


1.105820 


3.500000 


0.000000 
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Table 8: Control points for the optimal shape: fixed number of degrees of freedom case and cubic approximation is 
used in the B-splines representation. Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


3 


X 


y 


z 


i 


J 


X 


y 


z 








0.181755 


0.000000 


0.000000 


3 





-0.261467 


2.000000 


0.000000 





1 


0.318130 


0.000000 


0.089250 


3 


1 


-0.099029 


2.000000 


0.089250 





2 


0.454504 


0.000000 


0.126467 


3 


2 


0.063409 


2.000000 


0.126467 





3 


0.590878 


0.000000 


0.034380 


3 


3 


0.225847 


2.000000 


0.034380 





4 


0.727252 


0.000000 


0.091453 


3 


4 


0.388285 


2.000000 


0.091453 





5 


0.863626 


0.000000 


0.031200 


3 


5 


0.550723 


2.000000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.713161 


2.000000 


0.000000 


1 





0.240014 


0.333333 


0.000000 


4 





-0.009370 


2.666670 


0.000000 


1 


1 


0.359515 


0.333333 


0.089250 


4 


1 


0.148156 


2.666670 


0.089250 


1 


2 


0.479017 


0.333333 


0.126467 


4 


2 


0.305682 


2.666670 


0.126467 


1 


3 


0.598518 


0.333333 


0.034380 


4 


3 


0.463208 


2.666670 


0.034380 


1 


4 


0.718020 


0.333333 


0.091453 


4 


4 


0.620734 


2.666670 


0.091453 


1 


5 


0.837521 


0.333333 


0.031200 


4 


5 


0.778260 


2.666670 


0.031200 


1 


6 


0.957023 


0.333333 


0.000000 


4 


6 


0.935786 


2.666670 


0.000000 


2 





0.138733 


1.000000 


0.000000 


5 





0.198005 


3.500000 


0.000000 


2 


1 


0.250931 


1.000000 


0.089250 


5 


1 


0.353562 


3.500000 


0.089250 


2 


2 


0.363129 


1.000000 


0.126467 


5 


2 


0.509119 


3.500000 


0.126467 


2 


3 


0.475326 


1.000000 


0.034380 


5 


3 


0.664675 


3.500000 


0.034380 


2 


4 


0.587524 


1.000000 


0.091453 


5 


4 


0.820232 


3.500000 


0.091453 


2 


5 


0.699721 


1.000000 


0.031200 


5 


5 


0.975788 


3.500000 


0.031200 


2 


6 


0.811919 


1.000000 


0.000000 


5 


6 


1.131340 


3.500000 


0.000000 
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show mild changes when considering higher-degree polynomials and converge at the third degree. 
The optimal shape yields an increase in the lift while the overall wing area was kept the same as 
the baseline case. The required input power increases as well however. This indicates that not only 
changing the wing shape allows to achieve an efficiency of T] = 0.441 (~ 100% increase) but also 
more power needs to be introduced to the flying system. The set of optimal shapes obtained for 
the present case enables generation of higher thrust and lower lift but requires more aerodynamic 
power while maintaining almost the same level of efficiency in comparison with the single knot 
span case. This observation can be related to the wing area distribution. In fact, the proportion 
of area in the outer part of the wing (near the tip) which undergoes higher acceleration is thicker 
and then it needs more power to move. Besides, we remark that the optimal wings are longer 
(i.e., higher aspect ratio). The optimal wing tip, which is subjected to higher angular accelera- 
tion, is pushed further away. This explains the significant increase in the aerodynamic power with 
respect the baseline case. These observations provide insights on how to design efficient wing 
planforms based on the mission requirements. As such, shortening the wing tip helps in increas- 
ing the lift force and decreasing the required aerodynamic power. This feature can be useful for 
take-off flights. Furthermore, the differences observed in the shapes obtained from the different 
optimization cases explains the diversity of wing morphologies observed in nature. 

Table 9: Baseline vs. optimal results obtained for the fixed number of degrees of freedom case (k = 0.1). Lift, thrust, 
and area constraints are imposed. 



Wing shape 




V 


77 




p* 


A/A hl 


Baseline shape 




0.219 


4.484 


0.185 


0.844 


1 


Optimal shapes 


Linear 


0.439 


5.390 


0.498 


1.134 


1 




Quadratic 


0.440 


5.374 


0.506 


1.148 


1 




Cubic 


0.441 


5.398 


0.506 


1.146 


1 




Quartic 


0.441 


5.398 


0.506 


1.146 


1 



Figure [9] shows the lift, thrust, and aerodynamic power that develop over the flapping wing at 
/c — 0.1 for both the baseline and optimal cases, plotted as function of the flapping angle 0. As 
expected, the bulk of the useful aerodynamic forces (positive lift and thrust) are generated during 
the downstroke. In particular, thrust and lift peaks are reached near the middle of the downstroke 
phase. Positive lift is produced during both strokes, as the angle of attack induced by the flapping 
motion is smaller than the fixed pitch angle (5°) at the wing root, and thus the lift remains positive 
through almost the entire cycle. The optimal shape does not alter the phases of the aerodynamic 
quantities, since the effective angle of attack does not vary significantly over the flapping cycle. 
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Figure 9: Lift, thrust, and power plotted versus the flapping angle (j) for the baseline (rectangular) and optimal shapes 
= 0.1). 

Nevertheless, it is able to increase the time-averaged lift, thrust, and power (as seen in Table [9]), 
and also their peaks as shown in Figure [9] 

The vorticity in the wake was generated on and shed from the wing at an earlier time. As 
such, the wake is usually referred as the "historian" of the flow. Thus, examining the wake pattern 
and vorticity distribution is helpful to gain insight into the reasons why the obtained optimized 
shapes produce efficient flapping flights. The vorticity circulation strength of the wakes obtained 



for the baseline and optimal wing shapes is shown in Figure 10 Clearly, the overall strength of the 
wake has increased in comparison with the baseline case, as the average aerodynamic power has 
increased as shown in Table|9j In particular, stronger pockets of high circulation are observed in the 
wake aft of the optimal shape at the middle of the down and upstrokes (i.e., 4> ~ 0°). Furthermore, 
the vortex tip swirl is more pronounced for the baseline shape. So, the optimal shape managed to 
reduce the tip vortex effect which then produces higher thrust. 
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(a) Baseline case 



(b) Optimal case 



Figure 10: Wake patterns of the flapping wing obtained for the baseline and optimal shapes. Contour color levels 
denote the vorticity strength of the wing and the wake. The wake has been moved back for the sake of differentiating 
it from the wing. 



The unsteady pressure difference between the lower and upper surfaces of the wing for the 



baseline and optimal cases is plotted in Figure 1 1 The pressure difference varies over the flapping 



cycle, giving rise to a positive pressure jump during the downstroke and negative pressure jump 
during the upstroke. Clearly, the highest pressure levels are obtained during the downstroke. This 
explains the generation of the aerodynamic loads along the two phases of the flapping cycle as 
shown in Figure [9] Furthermore, higher pressure levels are achieved for the optimal case as reflec- 
tive of the increase in the required aerodynamic power. For both cases (baseline and optimal), we 
observe a pressure peak along the leading edge. This pressure peak may lead to higher adverse 
pressure gradient and may be indicative of an eventual flow separation ll30l . This peak is mostly 
related to the fact that the aerodynamic model (based on potential flow) enforces the flow to be 
attached around the sharp leading edge of the wing. 

Separate leading and trailing edges curvature effects. The above results indicate that one needs 
to design properly the trailing edge also since the wing may benefit from the interaction with the 
near wake and then improve its aerodynamic performance. In an attempt to separate the individual 
effects of the trailing and leading edges on the wing performance, we kept one of them to be 
fixed at 90° angle with respect to the flight direction while allowing the other to be curved and 
turned on the optimizer. Again, the number of degrees of freedom is held constant at 12 and the 
cubic approximation is used in the B-spline representation. Only increases of 17% and 1% in the 
propulsive efficiency were obtained for the optimal shapes when curving the leading and trailing 
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Pressure 



Flight direction 



-2.5 

(a) Baseline case 



»lml(«l 



Pressure 

^^^^^^^^^^^^^ Flight direction 



-2.5 5 
(b) Optimal case 

Figure 11: Pressure difference distribution over the wing surface during a flapping cycle. Top view of the flapping 
flight is shown. 

edges, separately. Clearly, both of them need to be varied simultaneously to enable superior design 
performances. 

Power constraint effect. The optimal shapes enable a significant improvement in terms of propul- 
sive efficiency, nevertheless, they require much more aerodynamic power to perform their forward 
flight. Thus, to keep the same level of aerodynamic power as required for the baseline case, a 
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power constraint is added and the optimization is reformulated as follows: 



Maximize i](x) 
Subject to L*(x) > L* hh 



A(x) < A hh (4) 



P*(x) < P* hh 

max \9i Ax) — f | < 9 CI , 

x e x, 



Results for the efficiency r\, lift L*, thrust T*, aerodynamic power P*, and area ratio A/A h \ 



obtained from the optimization case when imposing a power constraint are presented in Table 10 
A lower thrust has been obtained in comparison with the power unconstrained optimization case 
while maintaining almost the same levels of efficiency and lift. The decrease in the power con- 
sumption with respect to previous optimization case came at the price of a loss in the thrust force. 



Table 10: Baseline vs. optimal results, k = 0.1. Lift, thrust, area, and power constraints are imposed. Cubic 
approximation is used in the B-spline representation. 



Wing shape 


V 


L* 




p* 


A/A u 


Baseline shape 
Optimal shape 


0.219 
0.438 


4.484 
5.470 


0.185 
0.369 


0.844 
0.844 


1 

0.82 



To satisfy the power constraint, less area is distributed to the outer part of the optimal wing 



as shown in Figure 12 while the geometric parameters defining the optimal shapes are listed in 
Table \TT\ In particular, the wing is significantly shrunk near the tip. Since the acceleration resulting 
from the flapping motion is proportional to the distance from the wing root, the outer part of the 
wing (closer to the tip) undergoes greater accelerations and then contributes relatively more to the 
aerodynamic power requirement. Thus, to enable reduction in power consumption (with respect 
to the previous optimal configuration) while keeping the flapping flight as efficient as possible, the 
proportion of the wing near the tip becomes narrower. The wing's aspect ratio has increased in this 
power constrained optimal configuration as it did in the previous case, where the limiting factors 
were the wing area and the allowable range of variation adopted for the control points. 
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Table 11: Control points for the optimal shape: fixed number of degrees of freedom case and cubic approximation is 
used in the B-splines representation (k = 0.1). Lift, thrust, area, and power constraints are imposed. 

Control Points Control Points 



l 


J 


X 


V 


z 


i 


J 


X 


y 


z 








0.085173 


0.000000 


0.000000 


3 





-0.091147 


2.000000 


0.000000 





1 


0.237644 


0.000000 


0.089250 


3 


1 


0.031862 


2.000000 


0.089250 





2 


0.390115 


0.000000 


0.126467 


3 


2 


0.154870 


2.000000 


0.126467 





3 


0.542586 


0.000000 


0.034380 


3 


3 


0.277879 


2.000000 


0.034380 





4 


0.695058 


0.000000 


0.091453 


3 


4 


0.400887 


2.000000 


0.091453 





5 


0.847529 


0.000000 


0.031200 


3 


5 


0.523895 


2.000000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.646904 


2.000000 


0.000000 


1 





0.006480 


0.333333 


0.000000 


4 





0.183726 


2.666670 


0.000000 


1 


1 


0.165743 


0.333333 


0.089250 


4 


1 


0.254317 


2.666670 


0.089250 


1 


2 


0.325006 


0.333333 


0.126467 


4 


2 


0.324909 


2.666670 


0.126467 


1 


3 


0.484269 


0.333333 


0.034380 


4 


3 


0.395501 


2.666670 


0.034380 


1 


4 


0.643532 


0.333333 


0.091453 


4 


4 


0.466093 


2.666670 


0.091453 


1 


5 


0.802795 


0.333333 


0.031200 


4 


5 


0.536685 


2.666670 


0.031200 


1 


6 


0.962058 


0.333333 


0.000000 


4 


6 


0.607276 


2.666670 


0.000000 


2 





-0.233239 


1.000000 


0.000000 


5 





0.383476 


3.500000 


0.000000 


2 


1 


-0.042933 


1.000000 


0.089250 


5 


1 


0.460540 


3.500000 


0.089250 


2 


2 


0.147373 


1.000000 


0.126467 


5 


2 


0.537604 


3.500000 


0.126467 


2 


3 


0.337678 


1.000000 


0.034380 


5 


3 


0.614668 


3.500000 


0.034380 


2 


4 


0.527984 


1.000000 


0.091453 


5 


4 


0.691732 


3.500000 


0.091453 


2 


5 


0.718290 


1.000000 


0.031200 


5 


5 


0.768796 


3.500000 


0.031200 


2 


6 


0.908595 


1.000000 


0.000000 


5 


6 


0.845860 


3.500000 


0.000000 
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Flight direction 



Baseline No power constraint Power constraint 



Figure 12: Optimal wing shapes (with power constraint). Coordinates of the control points parameterizing the optimal 
shapes are reported in Table [TT] 

Aspect ratio effect. To gain a better understanding of the reasons for which the obtained optimized 
shapes enhance the performance of the flapping flight, the impact of varying the aspect ratio is 



examined for a rectangular wing. Figure 13 shows the sensitivity of the efficiency and power 



coefficient, defined as Cp = ?} T3 , to the wing's aspect ratio over a wide range of variation. We 
observe that the efficiency increases significantly with the aspect ratio to achieve a maximum near a 
value of 20, while the aerodynamic power increases monotonically when the wing's aspect ratio is 
increased. These results show that high- aspect-ratio wings generate more thrust and are more effi- 
cient than low-aspect-ratio wings however their flapping motion requires more aerodynamic power. 
This explains the identification of longer wings in the optimized configurations (the increase in the 
aspect ratio is limited by the area constraint and the allowable range of variation adopted for the 
control points) and because of the power constraint, the proportion of area in the outer region the 
wing has been reduced. This outer area is subjected to higher angular accelerations. 

Wings with higher aspect ratio enable better aerodynamic performance (when flying at mod- 
erate frequencies), as can be concluded from the present optimization study, however, such wings 
undergo greater deflections. This requires stricter structural design specifications to avoid unde- 
sirable aeroelastic effects, such as occurrence of flutter at low flight speeds. Besides, long wings 
have higher moment of inertia to overcome and create more parasitic drag ll5Dl which degrades the 
flight performance. All of the aforementioned reasons, which are not accounted for in our analysis, 
limit the consideration of high-aspect ratio wings when designing flying vehicles. 
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Figure 13: Variations of the efficiency and power coefficient with the wing's aspect ratio (AR) for a rectangular wing 
flapping at a reduced frequency k = 0.1. 



4.4.3. High reduced-frequency case 

We consider a higher reduced frequency (k = 0.4). Again, an aerodynamic mesh of 24 chord- 
wise elements and 20 spanwise elements and 120 time steps per flapping cycle are used. This 
configuration is found appropriate to achieve UVLM solution independency of wing and time dis- 
cretization. The previous case (« = 0.1) represents a quasi-steady flight condition, whereas in the 
current case unsteady effects are much more pronounced. The objective is to examine the impact 
of the reduced frequency on the optimal shapes. As expected, flying at higher reduced frequency 
involves more aerodynamic power. In fact, increasing the reduced frequency k leads to higher an- 
gular acceleration for the wing and consequently for the surrounding fluid (noncirculatory), which 
would require additional power to accelerate it. This noncirculatory fluid is pushed by the wing 
and leads to what is referred as the added mass effects. Looking at the optimal shapes obtained for 
the two reduced frequencies given in Figure [T4l we observe that the high-frequency wing shape is 



more tilted close the wing tip, resulting in a net drop in the area (as seen in Table 13). The geo- 



metric parameters controlling the optimal shapes are listed in Table 12 The optimal values of the 



lift, thrust, aerodynamic power, and area ratio are reported in Table 13 A moderate improvement 
in terms of propulsive efficiency (~ 28% increase) is obtained. Optimal results show different 
trends with respect to the low-frequency case. In fact, the area constraint is not active anymore. 
The thrust constraint is active and a lower aerodynamic power are obtained such that propulsive 
efficiency improves (from 0.457 to 0.576). Unlike the low-frequency case, the improvement in the 
propulsive efficiency is entirely due to a decrease in the input aerodynamic power. 



The lift, thrust, and power phase plots obtained at k = 0.4 are shown in Figure 15 Again, the 
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Table 12: Control points for the optimal shape: fixed number of degrees of freedom case and cubic approximation is 
used in the B-splines representation (k = 0.4). Lift, thrust, and area constraints are imposed. 

Control Points Control Points 



l 


3 


X 


V 


z 


i 


J 


X 


y 


z 








0.000227 


0.000000 


0.000000 


3 





0.030328 


2.000000 


0.000000 





1 


0.166856 


0.000000 


0.089250 


3 


1 


0.120033 


2.000000 


0.089250 





2 


0.333484 


0.000000 


0.126467 


3 


2 


0.209738 


2.000000 


0.126467 





3 


0.500113 


0.000000 


0.034380 


3 


3 


0.299443 


2.000000 


0.034380 





4 


0.666742 


0.000000 


0.091453 


3 


4 


0.389149 


2.000000 


0.091453 





5 


0.833371 


0.000000 


0.031200 


3 


5 


0.478854 


2.000000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.568559 


2.000000 


0.000000 







0.107023 


0.333333 


0.000000 


4 





0.256388 


2.666670 


0.000000 




1 


0.258868 


0.333333 


0.089250 


4 


1 


0.329331 


2.666670 


0.089250 




2 


0.410713 


0.333333 


0.126467 


4 


2 


0.402274 


2.666670 


0.126467 




3 


0.562558 


0.333333 


0.034380 


4 


3 


0.475218 


2.666670 


0.034380 




4 


0.714404 


0.333333 


0.091453 


4 


4 


0.548161 


2.666670 


0.091453 




5 


0.866249 


0.333333 


0.031200 


4 


5 


0.621104 


2.666670 


0.031200 




6 


1.018090 


0.333333 


0.000000 


4 


6 


0.694047 


2.666670 


0.000000 


2 





0.145300 


1.000000 


0.000000 


5 





0.433250 


3.398740 


0.000000 


2 


1 


0.262970 


1.000000 


0.089250 


5 


1 


0.510371 


3.398740 


0.089250 


2 


2 


0.380640 


1.000000 


0.126467 


5 


2 


0.587493 


3.398740 


0.126467 


2 


3 


0.498310 


1.000000 


0.034380 


5 


3 


0.664615 


3.398740 


0.034380 


2 


4 


0.615980 


1.000000 


0.091453 


5 


4 


0.741737 


3.398740 


0.091453 


2 


5 


0.733650 


1.000000 


0.031200 


5 


5 


0.818858 


3.398740 


0.031200 


2 


6 


0.851320 


1.000000 


0.000000 


5 


6 


0.895980 


3.398740 


0.000000 
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Flight direction 




Baseline Low frequency High frequency 



Figure 14: Effect of reduced frequency on optimal wing shapes. Coordinates of the control points parameterizing the 
optimal shapes are reported in Table [T2| 

Table 13: Baseline vs. optimal results, k = 0.4. Lift, thrust, and area constraints are imposed. Cubic approximation is 
used in the B-spline representation. 



Wing shape 


V 


L* 




p* 


A/Ah 


Baseline shape 
Optimal shape 


0.457 
0.576 


6.227 
7.294 


4.108 
4.108 


8.98 
7.13 


1 

0.716 



majority of the positive lift is generated during the downstroke, however now it reaches negative 
values during the upstroke. This is caused by the fact that the contribution of the flapping motion 
to the effective angle of attack is larger than that of the fixed pitch angle at the wing root, resulting 
in negative angles during the upstroke. 

The wakes trailing the baseline and optimal wing configurations obtained for high frequency 
(k = 0.4) are shown in Figure [To] where the circulation strength of those wakes can be observed. 
The vortex rings in the wake developed over a flapping cycle are located very close to the wing. 
Then, the wing-wake interactions have more significant effects on the aerodynamic behavior of the 
flapping wing. In comparison with the baseline case, the vortex rings in the wake aft of the optimal 
flapping wing are weak as reflective of the drop in required input power. 

4.5. Camber effect 

We allow the camber line to vary by introducing perturbations to the interior control points 
shown in Figure [JJ along the z-direction and turn on the optimizer. The camber line is kept the 
same at all cross sections along the spanwise direction. Moderate variations are considered to avoid 
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Figure 15: Lift, thrust, and power plotted versus the flapping angle <j> for the baseline (rectangular) and optimal shapes 
(/c = 0.4). 
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(a) Baseline case 



(b) Optimal case 



Figure 16: Wake patterns of the flapping wing obtained for the baseline and optimal shapes (k = 0.4). Contour 
color levels denote the vorticity strength of the wing and the wake. The wake has been moved back for the sake of 
differentiating it from the wing. 



wing cross sections with high angle of attack (with respect to the incoming freestream) that may 
initiate flow separation which cannot be captured in the present aerodynamic model. Again, the 
use of B-spline representation enables to smoothly discretize wing camber line with few degrees of 
freedom and then a variety of camber lines can be obtained rather than sticking to NACA profiles. 
The optimal results for the efficiency and aerodynamic loads and power are summarized in Tables 
14] and [15] for k — 0.1 and k = 0.4, respectively. As expected, the camber greatly affects the 



aerodynamic characteristics such as power and thrust. At low frequency, significantly higher thrust 
and efficiency can be obtained by optimizing the camber shape. The average value of the lift force 
in the optimized shape, obtained when introducing variations to the camber line, remained equal 
to the value in the baseline case. In fact, the optimizer would decrease the lift if it could, because 
this constraint is always active. On the other hand, higher thrust is achieved in the optimized 
configuration. So, there is a natural trade-off between lift and efficiency, but not so for thrust and 
efficiency. At high frequency, the optimization drives the design to a shape configuration with 
slightly lower input aerodynamic power, though the time-averaged thrust remained equal to that 
obtained for the baseline case, yielding a mild raise in the propulsive efficiency in comparison with 
the previous optimal shape obtained without varying the camber line. 

The optimal wing shapes and their camber lines are shown in Figures [FT] and [T8] for k = 0.1 



and k = 0.4, respectively, while the geometric configurations are detailed in Tables 16 and 17 The 
optimal shapes show qualitative similarity along the spanwise direction with respect to the previous 
optimal shape obtained without changing the camber line. The optimal camber line is thinner, in 
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particular for the low frequency case, than the NACA 83XX profile resulting in lower effective 
angle of attack. The change in the camber line enabled better performance in terms of efficiency 
but induced a reduction in lift generation (although the lift is higher than its value obtained for 
a rectangular wing). These results provide guidance on how to select the wing planform and the 
camber based on the nature of the mission and its aerodynamic requirements. 

Table 14: Baseline vs. optimal results, k = 0.1. Lift, thrust, and area constraints are imposed. Cubic approximation is 
used in the B-spline representation. 



Wing shape 


V 


L* 




p* 


A/A u 


Baseline shape 


0.219 


4.484 


0.185 


0.844 


1 


Optimal shape without camber variations 


0.441 


5.398 


0.506 


1.146 


1 


Optimal shape with camber variations 


0.533 


4.484 


0.526 


0.986 


0.93 



Table 15: Baseline vs. optimal results, k = 0.4. Lift, thrust, and area constraints are imposed. Cubic approximation is 
used in the B-spline representation. 



Wing shape 


V 


L* 




p* 


A/Au 


Baseline shape 


0.457 


6.227 


4.108 


8.98 


1 


Optimal shape without camber variations 


0.576 


7.294 


4.108 


7.13 


0.716 


Optimal shape with camber variations 


0.612 


6.910 


4.108 


6.499 


0.693 



Our study gives comparable optimal efficiencies to those obtained from previous studies HI 
[61 where the improvement of flight performance was achieved by introducing a dynamic shape 
deformation (usually referred as active shape morphing) which is based on imposing a prescribed 
bending and twisting. This shows again the importance of wing planform selection for an efficient 
design of flapping wing vehicles. Nevertheless, the optimal shapes reported here do not need to be 
morphed during forward flight to achieve their improved performance. 

5. Summary and Conclusion 

The wing shape/planform presents a critical determinant of performance in flapping flight. In 
this work, we use a three-dimensional version of the unsteady vortex lattice method (UVLM) to 
simulate the aerodynamic response of a flapping wing in forward flight and adopt B-spline rep- 
resentation to model the wing geometry. First, we carry out a convergence analysis to determine 
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Table 16: Control points for the optimal shape: fixed number of degrees of freedom case and cubic approximation 
is used in the B-splines representation (k = 0.1 and with camber variations). Lift, thrust, and area constraints are 
imposed. 

Control Points Control Points 



l 


J 


X 


y 


z 


i 


J 


X 


y 


z 








0.003370 


0.000000 


0.000000 


3 





-0.199021 


2.000000 


0.000000 





1 


0.169475 


0.000000 


0.089250 


3 


1 


-0.059942 


2.000000 


0.089250 





2 


0.335580 


0.000000 


0.042377 


3 


2 


0.079137 


2.000000 


0.042377 





3 


0.501685 


0.000000 


0.035764 


3 


3 


0.218216 


2.000000 


0.035764 





4 


0.667790 


0.000000 


0.053259 


3 


4 


0.357295 


2.000000 


0.053259 





5 


0.833895 


0.000000 


0.031200 


3 


5 


0.496374 


2.000000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.635453 


2.000000 


0.000000 


1 





0.091649 


0.333333 


0.000000 


4 





0.044826 


2.666670 


0.000000 


1 


1 


0.236577 


0.333333 


0.089250 


4 


1 


0.170684 


2.666670 


0.089250 


1 


2 


0.381505 


0.333333 


0.042377 


4 


2 


0.296542 


2.666670 


0.042377 


1 


3 


0.526433 


0.333333 


0.035764 


4 


3 


0.422401 


2.666670 


0.035764 


1 


4 


0.671361 


0.333333 


0.053259 


4 


4 


0.548259 


2.666670 


0.053259 


1 


5 


0.816289 


0.333333 


0.031200 


4 


5 


0.674117 


2.666670 


0.031200 


1 


6 


0.961217 


0.333333 


0.000000 


4 


6 


0.799975 


2.666670 


0.000000 


2 





0.001093 


1.000000 


0.000000 


5 





0.262770 


3.500000 


0.000000 


2 


1 


0.178065 


1.000000 


0.089250 


5 


1 


0.376402 


3.500000 


0.089250 


2 


2 


0.355037 


1.000000 


0.042377 


5 


2 


0.490034 


3.500000 


0.042377 


2 


3 


0.532009 


1.000000 


0.035764 


5 


3 


0.603666 


3.500000 


0.035764 


2 


4 


0.708981 


1.000000 


0.053259 


5 


4 


0.717298 


3.500000 


0.053259 


2 


5 


0.885953 


1.000000 


0.031200 


5 


5 


0.830930 


3.500000 


0.031200 


2 


6 


1.062930 


1.000000 


0.000000 


5 


6 


0.944562 


3.500000 


0.000000 
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Table 17: Control points for the optimal shape: fixed number of degrees of freedom case and cubic approximation 
is used in the B-splines representation (k = 0.4 and with camber variations). Lift, thrust, and area constraints are 
imposed. 

Control Points Control Points 



l 


J 


X 


y 


z 


i 


3 


X 


y 


z 








0.028839 


0.000000 


0.000000 


3 





-0.043142 


2.000000 


0.000000 





1 


0.190699 


0.000000 


0.089250 


3 


1 


0.061513 


2.000000 


0.089250 





2 


0.352559 


0.000000 


0.115212 


3 


2 


0.166168 


2.000000 


0.115212 





3 


0.514420 


0.000000 


0.021558 


3 


3 


0.270823 


2.000000 


0.021558 





4 


0.676280 


0.000000 


0.072092 


3 


4 


0.375477 


2.000000 


0.072092 





5 


0.838140 


0.000000 


0.031200 


3 


5 


0.480132 


2.000000 


0.031200 





6 


1.000000 


0.000000 


0.000000 


3 


6 


0.584787 


2.000000 


0.000000 


1 





0.010549 


0.333333 


0.000000 


4 





0.285517 


2.666670 


0.000000 


1 


1 


0.183619 


0.333333 


0.089250 


4 


1 


0.337391 


2.666670 


0.089250 


1 


2 


0.356689 


0.333333 


0.115212 


4 


2 


0.389265 


2.666670 


0.115212 


1 


3 


0.529759 


0.333333 


0.021558 


4 


3 


0.441138 


2.666670 


0.021558 


1 


4 


0.702830 


0.333333 


0.072092 


4 


4 


0.493012 


2.666670 


0.072092 


1 


5 


0.875900 


0.333333 


0.031200 


4 


5 


0.544886 


2.666670 


0.031200 


1 


6 


1.048970 


0.333333 


0.000000 


4 


6 


0.596760 


2.666670 


0.000000 


2 





0.401749 


1.000000 


0.000000 


5 





0.474763 


3.500000 


0.000000 


2 


1 


0.482217 


1.000000 


0.089250 


5 


1 


0.537329 


3.500000 


0.089250 


2 


2 


0.562685 


1.000000 


0.115212 


5 


2 


0.599895 


3.500000 


0.115212 


2 


3 


0.643152 


1.000000 


0.021558 


5 


3 


0.662461 


3.500000 


0.021558 


2 


4 


0.723620 


1.000000 


0.072092 


5 


4 


0.725027 


3.500000 


0.072092 


2 


5 


0.804088 


1.000000 


0.031200 


5 


5 


0.787593 


3.500000 


0.031200 


2 


6 


0.884556 


1.000000 


0.000000 


5 


6 


0.850159 


3.500000 


0.000000 
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Figure 17: Optimal wing shape and camber line (k = 0.1). Cubic approximation is used in the B-spline representation. 



Coordinates of the control points parameterizing the optimal shapes are reported in Table 16 



the appropriate time step and aerodynamic mesh sizes which are needed to reach UVLM solu- 
tion independence to discretization. In particular, selecting the proper time step was important 
to observe the UVLM solution when refining the aerodynamic mesh. Then, we combine UVLM 
with a gradient-based optimizer to identify a set of efficient shapes that optimizes the aerodynamic 
performance. Flow simulations using the optimal wing shapes indicate that changes in the shape 
have significant effects on relevant averaged quantities such as lift, thrust, and aerodynamic power. 
Increasing the number of variables (i.e., providing the wing shape with a greater degree of spatial 
freedom) enables increasingly superior designs. The optimization study shows that the camber line 
and the leading and trailing edges represent the key wing shape parameters that control the gen- 
eration of the aerodynamic loads and flight performance depending on the nature of the mission 
assigned to the flying vehicle. Furthermore, the optimal shapes show significant dependency on 
the reduced frequency. At low frequency, the optimization pushes the design to a shape configura- 
tion with substantial increase in the time-averaged thrust, while the average aerodynamic power is 
increased, resulting in a significant increase in the propulsive efficiency. To keep the same level of 
aerodynamic power as required for the baseline rectangular wing, a power constraint is added to 
the optimization problem. A reduction in the outer region of the wing near the tip is observed in 
the optimal shape. This outer area is subjected to higher angular accelerations. At high frequency, 
the optimization drives the design to a shape configuration with lower input aerodynamic power, 
though the time-averaged thrust remained equal to that obtained for the baseline case, resulting in 
a moderate improvement in terms of propulsive efficiency. 

This work was concerned solely with the assessment of the aerodynamic performance of rigid 
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Flight direction 




Camber line: baseline vs. optimal 

Baseline Optimal Optimal 

(without camber (with camber 
variations) variations) 



Figure 18: Optimal wing shape and camber line (n = 0.4). Cubic approximation is used in the B-spline representation. 



Coordinates of the control points parameterizing the optimal shapes are reported in Table 17 



flapping wings when considering different shapes. Implementing a full aeroelastic framework that 
couples a nonlinear shell model and UVLM to test the obtained optimal shapes and check their 
superiority over the baseline configuration is the topic of our current research effort. 
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